Für die erste Abgabe im Modul PRR beschäftigen wir uns mit (multiplen) linearen und logistischen Regressionen. Diese Funktionen wenden wir auf den Datensatz “Carseats” an. Wir haben uns für diese Methoden entschieden, weil sie ein zentraler Bestandteil von Data Science sind. Sie wird verwendet, um Zusammenhängen in Datensätze zu verstehen und um Prognosen erstellen zu können. Diese Regressionen sind im R Grundpakte enthalten. Es muss also keine separate Library installiert und geladen werden.
Carseats ist ein generischer Datensatz, der im Buch “An Introduction to Statistical Learning” verwendet wird. Dieses Buch erarbeiten wir in diesem Studiengang im Modul “Statistisches Lernen”. Dort wird er zwar hauptsächlich angewendet, um Decision Trees zu erstellen. Dank 11 Variablen und 400 Observationen eignet er sich aber auch gut für Regressionen. Jede Zeile im Datensatz entspricht einem Verkaufslokal, gejoint mit Durchschnittsstatitistiken für das Einzugsgebiet des Lokals. Für eine genau Erklärung der Variablen verweisen wir auf die Homepage https://rdrr.io/cran/ISLR2/man/Carseats.html. Für die unten stehenden Hypothesen werden wir mit folgenden Variablen arbeiten:
Bei der Linearen Regression versucht man eine Zielvariable in Abhängigkeit durch eine andere Variable vorherzusagen. Das daraus entstehende Modell wird als eine Gerade dargestellt mit der
Formel: y = f(x) = ax + b = \(\beta_{0}\) + \(\beta_{1}\)x = mx + n
Das Wissen wird immersiv durch die Module llr, mlr und stl abgedeckt.
Bei der Linearen Multiple Regression versucht man eine Zielvariable in Abhängigkeit durch mehrere andere Variable vorherzusagen.
Formel ausgeschrieben: y = f(X) = \(\beta_{0}\) + \(\beta_{1}\) \(x_{1}\) + \(\beta_{2}\) \(x_{2}\) + … + \(\beta_{n}\) \(x_{n}\)
Formel vektorisiert: y = f(X) = X \(\cdot\) \(\vec{\beta_{n}}\) + \(\beta_{0}\)
Das Wissen wird immersiv durch die Module llr, mlr und stl abgedeckt.
Bei der Logistischen Regression versucht man eine binäre Zielvariable in Abhängigkeit durch eine andere Variable vorherzusagen.
Sigmoid Formel: \(\sigma(z) = \frac{1}{1 + e^{-z}}\)
Formel fuer logistische Regression: f(x) = \(\sigma(ax + b)\) = \(\frac{1}{1 + e^{-(ax + b)}}\)
Das Wissen wird immersiv durch die Module llr, mlr und stl abgedeckt.
Zur Analyse des Datensatzes haben wir uns folgende Fragen gestellt:
Sales ist abhängig von Income/Advertising/Population. (Lineare Regression)
Die Differenz von Comprice und Price hat einen Einfluss auf den Verkauf. (Lineare Regression)
Sales lässt sich durch die Verwendung von mehreren Features zuverlässiger vorhersagen. (Multiple Lineare Regression)
Bewohner von urbanen Regionen kaufen weniger Kindersitze als Bürger ländlicher Regionen.
In diesem Kapitel installieren wir die notwendigen Libraries und verschaffen uns einen Überlick über den Carseats Datensatz.
# Libraries installieren und laden
library(ISLR2)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0 ✔ purrr 0.3.5
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.1 ✔ stringr 1.5.0
## ✔ readr 2.1.3 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
In der Library ISLR2 befindet sich der Datensatz. Tidyverse wird für die Darstellung des Datensatzes und Visualisierungen benötigt.
In einem ersten Schritt ist es wichtig den Datensatz zu verstehen, bevor Visualisierungen und Modelle erstellt werden.
# Carseats Daten einlesen
data(Carseats)
# Die ersten 6 Oberservationen ausgeben
head(Carseats)
## Sales CompPrice Income Advertising Population Price ShelveLoc Age Education
## 1 9.50 138 73 11 276 120 Bad 42 17
## 2 11.22 111 48 16 260 83 Good 65 10
## 3 10.06 113 35 10 269 80 Medium 59 12
## 4 7.40 117 100 4 466 97 Medium 55 14
## 5 4.15 141 64 3 340 128 Bad 38 13
## 6 10.81 124 113 13 501 72 Bad 78 16
## Urban US
## 1 Yes Yes
## 2 Yes Yes
## 3 Yes Yes
## 4 Yes Yes
## 5 Yes No
## 6 No Yes
Carseats <- Carseats %>% select(Sales, Price, CompPrice, Advertising, ShelveLoc, Urban, US,
Population, Income, Age, Education,)
# Die letzten 6 Observationen ausgeben
tail(Carseats)
## Sales Price CompPrice Advertising ShelveLoc Urban US Population Income Age
## 395 5.35 139 130 19 Bad Yes Yes 366 58 33
## 396 12.57 128 138 17 Good Yes Yes 203 108 33
## 397 6.14 120 139 3 Medium No Yes 37 23 55
## 398 7.41 159 162 12 Medium Yes Yes 368 26 40
## 399 5.94 95 100 7 Bad Yes Yes 284 79 50
## 400 9.71 120 134 0 Good Yes Yes 27 37 49
## Education
## 395 16
## 396 14
## 397 11
## 398 18
## 399 12
## 400 16
Durch die Ausgabe von Kopfzeile und Fusszeile kommen wir mit dem Datensatz in Berührung und können uns ein Bild der Attribute sowie Observationen machen.
# Informationen zum DataFrame ausgeben
str(Carseats)
## 'data.frame': 400 obs. of 11 variables:
## $ Sales : num 9.5 11.22 10.06 7.4 4.15 ...
## $ Price : num 120 83 80 97 128 72 108 120 124 124 ...
## $ CompPrice : num 138 111 113 117 141 124 115 136 132 132 ...
## $ Advertising: num 11 16 10 4 3 13 0 15 0 0 ...
## $ ShelveLoc : Factor w/ 3 levels "Bad","Good","Medium": 1 2 3 3 1 1 3 2 3 3 ...
## $ Urban : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 2 2 1 1 ...
## $ US : Factor w/ 2 levels "No","Yes": 2 2 2 2 1 2 1 2 1 2 ...
## $ Population : num 276 260 269 466 340 501 45 425 108 131 ...
## $ Income : num 73 48 35 100 64 113 105 81 110 113 ...
## $ Age : num 42 65 59 55 38 78 71 67 76 76 ...
## $ Education : num 17 10 12 14 13 16 15 10 10 17 ...
Durch str wird uns die Grösse des Dataframes angezeigt, sowie die Anzahl an Features (Spalten). Weiter erkennen wir von jeder Spalte den Datenyp und auch einige Werte, die in der Spalte vorkommen.
Carseats$ShelveLoc <- factor(Carseats$ShelveLoc, ordered = TRUE, levels = c("Good", "Medium", "Bad"))
# Statistische Kennzahlen
summary(Carseats)
## Sales Price CompPrice Advertising ShelveLoc
## Min. : 0.000 Min. : 24.0 Min. : 77 Min. : 0.000 Good : 85
## 1st Qu.: 5.390 1st Qu.:100.0 1st Qu.:115 1st Qu.: 0.000 Medium:219
## Median : 7.490 Median :117.0 Median :125 Median : 5.000 Bad : 96
## Mean : 7.496 Mean :115.8 Mean :125 Mean : 6.635
## 3rd Qu.: 9.320 3rd Qu.:131.0 3rd Qu.:135 3rd Qu.:12.000
## Max. :16.270 Max. :191.0 Max. :175 Max. :29.000
## Urban US Population Income Age
## No :118 No :142 Min. : 10.0 Min. : 21.00 Min. :25.00
## Yes:282 Yes:258 1st Qu.:139.0 1st Qu.: 42.75 1st Qu.:39.75
## Median :272.0 Median : 69.00 Median :54.50
## Mean :264.8 Mean : 68.66 Mean :53.32
## 3rd Qu.:398.5 3rd Qu.: 91.00 3rd Qu.:66.00
## Max. :509.0 Max. :120.00 Max. :80.00
## Education
## Min. :10.0
## 1st Qu.:12.0
## Median :14.0
## Mean :13.9
## 3rd Qu.:16.0
## Max. :18.0
Wir geben von jeder numerischen Spalte die summarischen Statistiken aus. Darin enthalten sind: Minimum, 1. Quantil, Median, Mittelwert, 3. Quantil, Maximum. Im summary() werden auch die Anzahl fehlende Werte angezeigt. In unserem Fall besitzt der Datensatz keine fehlende Werte. Die Aufbereitung der Daten entfällt deshalb.
ggplot(data = Carseats, aes(x = Sales)) +
geom_density(fill = "lightgreen", color = "black", alpha = 0.3) +
labs(x = "Sales",
y = "Anzahl",
title = "Verteilung der Sales",
subtitle = "Carseats Datensatz")
In diesem Dichte-Diagramm erkennen wir die Verteilung der numerischen Variabel Sales. Die y Achse ist die Wahrscheinlichkeit für das Auftreten einer numerischen Sales Variabel. Aus diesem Grund wird die Dichtefunktion auch oft Wahrscheinlichkeitsdichtefunktion genannt. Sie ist eine gegelättete Version des Histogramm und wird oft nach dem gleichen Konzept verwendet.
ggplot(data = Carseats, aes(x = Sales, fill = Urban)) +
geom_density(alpha = 0.3) +
labs(x = "Sales",
y = "Anzahl",
title = "Verteilung der Sales",
subtitle = "Carseats Datensatz")
Analog zur ersten Visualisierung betrachten wir hier die
Dichteverteilung von Sales mit der Auftrennung von Urban (ja/nein),
sprich können wir mit diesem Plot unterschiede zwischen Urbanen Regionen
(städtische oder ländliche) Sales betrachten. In unserem Fall gibt es
keinen grosse Unterschied, man erkennt ganz leicht, dass Sales in nicht
ländlichen Regionen diese im Sales Bereich grösser ist als bei
Städtischen Regionen.
ggplot(Carseats, aes(x = Sales, fill = US)) +
geom_density(alpha = 0.3) +
facet_grid(Urban ~ ShelveLoc, labeller = label_both) +
labs(x = "Sales",
y = "Anzahl",
title = "Verteilung der Sales unterteilt nach ShelveLoc und Urban",
subtitle = "Carseats Datensatz")
In diesem Multiplot unterteilen wir die Dichteverteilung nach weiteren
Kategorien, wie Urban und ShelveLoc. Auch unterteilen wir farblich, ob
sich das Geschäft in den US befindet oder nicht, damit wollen wir ein
Gefühl für die Sales Verteilung mit verschiedenen Kategorischen Werten
kriegen. Aufällig ist, dass die Verteilungen bei ShelveLoc Bad deutlich
weiter Links sind (Mittelwert tiefer als bei ShelveLoc Medium und Good).
Auch erkennen wir bei ShelveLoc Good, dass die Kindersitze in den nicht
US Ländern weniger verkauft werden als in den US selber.
ggplot(Carseats, aes(x = Price, fill = US)) +
geom_density(alpha = 0.3) +
facet_grid(Urban ~ ShelveLoc, labeller = label_both) +
labs(x = "Price",
y = "Anzahl",
title = "Verteilung von Price unterteilt nach ShelveLoc und Urban",
subtitle = "Carseats Datensatz")
Analog zum Sales multiplot erstellen wir ein Multiplot von Price, um
eine Übersicht von Price mit unterschiedlichen Kategorien zu kriegen.
Wir erkennen, dass bei ShelveLoc Good und Urbanen Regionen die Preise in
den US deutlich höher sind als in nicht US Ländern. Auch erkennen wir
bei ShelveLoc Bad und nicht Urbanen Regionen, das die Preise für einen
Kindersitz bei nicht US Ländern teurer sind.
Thema: Einfache lineare Regression
plot(Sales ~ Income, data = Carseats,
main = "Sales vs Income", pch=21, bg = "lightblue")
Bevor wir mit der Linearen Regression beginnen, erstellen wir zuerst ein Streudiagramm von beiden Variabeln, um die Zusammenhänge zu visualiseren. Hier erkennen wir im Plot, dass die Variabeln Einkommen und Sales nicht besonders stark korrelieren, da die Punkte überall im Plot verstreut sind.
lm.fit1 <- lm(Sales ~ Income, data = Carseats)
Mit der Funktion “lm()” wird das Lineare Modell erstellt und trainiert.
summary(lm.fit1)
##
## Call:
## lm(formula = Sales ~ Income, data = Carseats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.2629 -1.9447 -0.1772 1.7654 8.9064
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.44356 0.37061 17.386 < 2e-16 ***
## Income 0.01533 0.00500 3.067 0.00231 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.795 on 398 degrees of freedom
## Multiple R-squared: 0.02309, Adjusted R-squared: 0.02063
## F-statistic: 9.407 on 1 and 398 DF, p-value: 0.00231
R hat nun die linear Regression erstellt. Die wichtigsten Werte lassen sich mit der Methode summary() anzeigen. Der Intercept beträgt 6.44 und die Steigung 0.015. Pro 1’000 Dollar höheres Einkommen steigt der Verkauf um 15 Stück. Dies erscheint im Verhältnis zum Intercept wenig, es müsste darüber diskutiert werden, ob diese Erkenntnis für einen möglichen Kunden aussagekräftig ist. Anhand vom P Wert (0.00231) lässt sich aber feststellen, dass die Steigung siginfikant ist. Das Einkommen der Region hat also einen Einfluss auf die Verkaufszahlen. Vor einer endgültigen Aussage müssen aber weitere Untersuchungen, wie die Überprüfung der Residuen, durchgeführt werden. Diese Untersuchung führen wir später durch.
names(lm.fit1)
## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "xlevels" "call" "terms" "model"
Das Modell enthält weitere Berechnungen und Kennzahlen. Names() zeigt eine Übersicht dieser Methoden.
lm.fit1$coefficients
## (Intercept) Income
## 6.44355748 0.01533361
Die oben ersichtlichen Koeffizienten lassen sich direkt via $coefficients abrufen.
confint(lm.fit1)
## 2.5 % 97.5 %
## (Intercept) 5.714962797 7.17215216
## Income 0.005504876 0.02516235
Wertvoll ist die Berechnung der Konfidenzintervalle, die im Summary nicht angezeigt wird. Diese gibt an, in welchem Bereich sich die Koeffizienten mit einer Wahrscheinlichkeit von 95% befinden.
plot(Sales ~ Income, data = Carseats,
main = "Sales vs Income", pch=21, bg = "lightblue")
abline(lm.fit1, col = "red")
Das Modell lässt sich direkt mit der Plotfunktion darstellen. Mit abline() lassen sich Linien zeichnen, in diesem Fall die lineare Regression. Im Plot erkennen wir, dass Income und Sales eine Punktwolke bilden. Sales streut auf der ganzen X Achse von maximalen bis minimalen beobachteten Wert (0 - 15). Entsprechend korrelieren beide Variabeln wenig miteinander. Es gibt keinen klaren Zusammenhang.
par(mfrow = c(2, 2))
plot(lm.fit1)
Nun, da wir unser Modell haben, können wir feststellen, ob es ein gutes Modell ist, indem wir bei den Residuen folgende 3 Bedingungen überprüfen:
Das mittlere Kriterium lässt sich aus der Summary Funktion ablesen. Der Median der Residuen entspricht -0.18, dieses Kriterium ist erfüllt.
Auch dafür bietet R einfache Möglichkeiten an. Mit der folgenden beiden Funktionen lassen sich vier unterschiedliche Überprüfungen der Residuen anzeigen.
Im ersten Plot oben links (Tukeyanscombe Plot) wird ersichtich, dass die Residuen/Fehler unabhängig von einander sind und nicht keinem Muster folgen. Auch dieses Kriterium ist somit erfüllt.
Die Verteilung der Residuen lässt sich einfach im Histogramm überprüfen, ob deren Mittelwert sich um 0 Verteilt.
Überprüfen vom Erwartungswert = 0 und Normalverteilt:
hist(lm.fit1$residuals,
main = "Verteilung der Residuen", col = "lightblue")
Hypothese: Sales ist abhängig von Income
Wir erkennen aufgrund des Streudiagramm, dass eine Korrelation zwischen Einkommen und Sales kaum vorhanden ist. Das Lineare Modell hat einen R^2 von 0.02, dies bestätigt die Aussage.
Aus diesem Grund ist Sales nicht abhängig von Income.
plot(Sales ~ Population, data = Carseats,
main = "Sales vs Population", pch=21, bg = "lightblue")
Bevor das Modell erstellt wird, visualisieren wir beide Variabeln Population und Sales in einem Streudiagramm.
lm.fit2 <- lm(Sales ~ Population, data = Carseats)
summary(lm.fit2)
##
## Call:
## lm(formula = Sales ~ Population, data = Carseats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.5864 -2.0176 -0.0597 1.6892 8.7213
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.2401837 0.2906658 24.909 <2e-16 ***
## Population 0.0009672 0.0009593 1.008 0.314
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.824 on 398 degrees of freedom
## Multiple R-squared: 0.002547, Adjusted R-squared: 4.116e-05
## F-statistic: 1.016 on 1 and 398 DF, p-value: 0.314
Die Einfluss der Population ist nicht signifikant. Der folgende Plot bestätigt, dass die Steigung 0 beträgt.
plot(Sales ~ Population, data = Carseats,
main = "Sales vs Population", pch=21, bg = "lightblue")
abline(lm.fit2, col = "red")
Hier in diesem Plot stellen wir die Population mit Sales sowie dem Modell zusammen. Wir erkennen, dass die Punkte sich deutlich von der roten Linie (Model) streuen. Das Problem besteht analog zu Sales und Income.
par(mfrow = c(2, 2))
plot(lm.fit2)
Auch in diesem Model überpruefen wir die folgenden Kriterien des Models.
Das mittlere Kriterium lässt sich aus der Summary Funktion ablesen. Der Median der Residuen entspricht -0.0597, dieses Kriterium ist erfüllt.
Auch dafür bietet R einfache Möglichkeiten an. Mit der folgenden beiden Funktionen lassen sich vier unterschiedliche Überprüfungen der Residuen anzeigen.
Im ersten Plot oben links (Tukeyanscombe Plot) wird ersichtich, dass die Residuen/Fehler unabhängig von einander sind und nicht einem Muster folgen. Auch dieses Kriterium ist somit erfüllt.
Die Verteilung der Residuen lässt sich einfach im Histogramm überprüfen und ob sich deren Mittelwert um 0 verteilt.
hist(lm.fit2$residuals,
main = "Verteilung der Residuen", col = "lightblue")
Hypothese: Sales ist abhängig von Population
Leider auch hier, ähnlich wie zu Sales vs. Income erkennen wir keinen Trend. Der R^2 Score beträgt 0.002547, es handelt sich also auch um kein verlässliches Modell.
Aus diesem Grund, stimmt auch diese Hypothese nicht.
plot(Sales ~ Advertising, data = Carseats,
main = "Sales vs Advertising", pch=21, bg = "lightblue")
Bevor das Modell erstellt wird, visualisieren wir die beiden Variabeln Advertising und Sales.
lm.fit3 <- lm(Sales ~ Advertising, data = Carseats)
summary(lm.fit3)
##
## Call:
## lm(formula = Sales ~ Advertising, data = Carseats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.3770 -1.9634 -0.1037 1.7222 8.3208
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.7370 0.1925 35.007 < 2e-16 ***
## Advertising 0.1144 0.0205 5.583 4.38e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.723 on 398 degrees of freedom
## Multiple R-squared: 0.07263, Adjusted R-squared: 0.0703
## F-statistic: 31.17 on 1 and 398 DF, p-value: 4.378e-08
plot(Sales ~ Advertising, data = Carseats,
main = "Sales vs Advertising", pch=21, bg = "lightblue")
abline(lm.fit3, col = "red")
Leider erkennen wir schon hier, dass es keinen Trend von Sales vs. Advertising vorhanden ist. Vollständigkeitshalber führen wir die Analysen trotzdem durch.
par(mfrow = c(2, 2))
plot(lm.fit3)
Auch in diesem Model überpruefen wir die folgenden Kriterien des Models.
Das mittlere Kriterium lässt sich aus der Summary Funktion ablesen. Der Median der Residuen entspricht -0.1037, dieses Kriterium ist erfüllt.
Auch dafür bietet R einfache Möglichkeiten an. Mit der folgenden beiden Funktionen lassen sich vier unterschiedliche Überprüfungen der Residuen anzeigen.
Im ersten Plot oben links (Tukeyanscombe Plot) wird ersichtich, dass die Residuen/Fehler unabhängig von einander sind und nicht keinem Muster folgen. Auch dieses Kriterium ist somit erfüllt.
Die Verteilung der Residuen lässt sich einfach im Histogramm überprüfen und ob sich deren Mittelwert um 0 verteilt.
hist(lm.fit3$residuals,
main = "Verteilung der Residuen", col = "lightblue")
Hypothese: Sales ist abhängig von Advertising
Leider ist auch bei Advertising und Sales kein linearer Zusammenhang erslichtlich. In der 3. Hypothese beschäftigen wir uns mit der multiplen Linearen Regression, wir erwarten, dass sie besser sein wird.
Keines der Linearen Modelle hat einen Linearen Zusammenhang. Mit einem einfachen Linearen Modell ist es nicht möglich, Sales vorherzusagen.
Thema: Einfach lineare Regression mit berechneten Daten
Hinzufügen einer Spalte mit der Differenz zwischen unserem Preis und der Konkurrenz.
Carseats <- Carseats %>% mutate(DiffPrice = Price - CompPrice)
str(Carseats)
## 'data.frame': 400 obs. of 12 variables:
## $ Sales : num 9.5 11.22 10.06 7.4 4.15 ...
## $ Price : num 120 83 80 97 128 72 108 120 124 124 ...
## $ CompPrice : num 138 111 113 117 141 124 115 136 132 132 ...
## $ Advertising: num 11 16 10 4 3 13 0 15 0 0 ...
## $ ShelveLoc : Ord.factor w/ 3 levels "Good"<"Medium"<..: 3 1 2 2 3 3 2 1 2 2 ...
## $ Urban : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 2 2 1 1 ...
## $ US : Factor w/ 2 levels "No","Yes": 2 2 2 2 1 2 1 2 1 2 ...
## $ Population : num 276 260 269 466 340 501 45 425 108 131 ...
## $ Income : num 73 48 35 100 64 113 105 81 110 113 ...
## $ Age : num 42 65 59 55 38 78 71 67 76 76 ...
## $ Education : num 17 10 12 14 13 16 15 10 10 17 ...
## $ DiffPrice : num -18 -28 -33 -20 -13 -52 -7 -16 -8 -8 ...
plot(Sales ~ DiffPrice, data = Carseats,
main = "Sales vs DiffPrice", pch=21, bg = "lightblue")
In dieser Visualisierung von DiffPrice und Sales erkennen wir, dass die Streupunkte einem gewissen Trend folgen.
lm.fit4 <- lm(Sales ~ DiffPrice, data = Carseats)
summary(lm.fit4)
##
## Call:
## lm(formula = Sales ~ DiffPrice, data = Carseats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5885 -1.6583 -0.2678 1.5394 6.2367
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.691590 0.125567 53.29 <2e-16 ***
## DiffPrice -0.087662 0.005891 -14.88 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.267 on 398 degrees of freedom
## Multiple R-squared: 0.3575, Adjusted R-squared: 0.3559
## F-statistic: 221.5 on 1 and 398 DF, p-value: < 2.2e-16
Die Auswertung zeigt, dass dieses Modell einen R^2 von 0.3575 hat. Dies ist nicht besonders gut, dennoch besser als die ersten drei Modelle.
plot(Sales ~ DiffPrice, data = Carseats,
main = "Sales vs DiffPrice", pch=21, bg = "lightblue")
abline(lm.fit4, col = "red")
Hier in diesem Diagramm erkennen wir, dass das Modell eine negative Steigung hat. Dies erkennen wir einerseits an der Steigung im Summary, aber auch bei der Visualisierung.
par(mfrow = c(2, 2))
plot(lm.fit4)
Auch in diesem Modell uüberpruefen wir die folgenden Kriterien des Modells.
Das mittlere Kriterium lässt sich aus der Summary Funktion ablesen. Der Median der Residuen entspricht -0.2678, dieses Kriterium ist erfüllt.
Auch dafür bietet R einfache Möglichkeiten an. Mit der folgenden beiden Funktionen lassen sich vier unterschiedliche Überprüfungen der Residuen anzeigen.
Im ersten Plot oben links (Tukeyanscombe Plot) wird ersichtich, dass die Residuen/Fehler unabhängig von einander sind und nicht keinem Muster folgen. Auch dieses Kriterium ist somit erfüllt.
Die Verteilung der Residuen lässt sich einfach im Histogramm überprüfen und ob sich deren Mittelwert um 0 verteilt.
hist(lm.fit4$residuals,
main = "Verteilung der Residuen", col = "lightblue")
Hypothese: Die Differenz von Comprice und Price hat einen Einfluss auf den Verkauf.
Die zweite Hypothese hat einen deutlich besseren R^2 Score verglichen zu den ersten drei Modellen bei der Hypothese 1. Die Differenz von Comprice und Price hat aber keinen grossen Einfluss auf Sales, weil das Bestimmtheitsmass 0.3575 beträgt. Weil sich die Residuen also nur teilweise durch das Modell erklären lassen, verwerfen wir auch die zweite Hypothese.
Es ist interessant zu erkennen, dass der Preisunterschied zwischen der Konkurrenz und dem Produkt einen ganz leichten negativen Trend zeigt. Nichts destotrotz, sind wir mit diesem Modell nicht zufrieden und begeben uns direkt zur Hypothese 3.
Thema: Multiple lineare Regression
Nachdem wir festgestellt haben, dass sich Sales nicht durch einzelne Variabeln erklären lasst, untersuchen wir die Möglichkeit, ob Sales sich durch mehrere Features erklären lässt. Wir werden sie zuerst gegenüber zwei Variabel untersuchen. Danach schauen wir uns den ganzen Datensatz an und behalten die aussagekräftigsten Variabeln.
In diesem Kapitel erstellen wir mehrere Modelle und wählen eines aufgrund einer ersten Auswertung für die weitere Untersuchung aus.
mlr.fit1 <- lm(Sales ~ Income + Population, data = Carseats)
summary(mlr.fit1)
##
## Call:
## lm(formula = Sales ~ Income + Population, data = Carseats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.0371 -1.9342 -0.1477 1.7786 8.8532
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.1785056 0.4493283 13.751 < 2e-16 ***
## Income 0.0153747 0.0049991 3.075 0.00225 **
## Population 0.0009902 0.0009493 1.043 0.29757
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.795 on 397 degrees of freedom
## Multiple R-squared: 0.02576, Adjusted R-squared: 0.02085
## F-statistic: 5.248 on 2 and 397 DF, p-value: 0.005627
Zuerst untersuchen wir Sales gegenüber Income und der Population. In der Auswertung wird ersichtlich, dass Population keinen Zusammenhang aufweist. R^2 spricht ebenfalls gegen einen Zusammenhang in dieser Auswertung.
Die F-Statistik sagt aus, dass mit hoher Wahrscheinlichkeit (p-Wert < 0.6%) mindestens eine Variabel einen Einfluss hat. In unserem Fall ist das Income, was wir schon aus vorhergehenden Modellen wissen.
mlr.fit2 <- lm(Sales ~ Income * Population, data = Carseats)
summary(mlr.fit2)
##
## Call:
## lm(formula = Sales ~ Income * Population, data = Carseats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.9097 -1.9657 -0.2132 1.7682 8.8761
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.176e+00 7.950e-01 9.027 <2e-16 ***
## Income 1.161e-03 1.060e-02 0.110 0.913
## Population -2.655e-03 2.578e-03 -1.030 0.304
## Income:Population 5.198e-05 3.419e-05 1.520 0.129
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.79 on 396 degrees of freedom
## Multiple R-squared: 0.03141, Adjusted R-squared: 0.02408
## F-statistic: 4.281 on 3 and 396 DF, p-value: 0.005445
Hier haben wir untersucht, ob die Variabeln Income und Population nicht nur additiv einen Einfluss haben, sondern ob sie auch gemeinsam wirken (Interaktion/Synergie). Das ist in diesem Fall aber nicht gegeben (P-Wert von 0.129 ist zu hoch).
mlr.fit3 <- lm(Sales ~ ., data = Carseats)
summary(mlr.fit3)
##
## Call:
## lm(formula = Sales ~ ., data = Carseats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8692 -0.6908 0.0211 0.6636 3.4115
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.9295889 0.5963763 13.296 < 2e-16 ***
## Price -0.0953579 0.0026711 -35.700 < 2e-16 ***
## CompPrice 0.0928153 0.0041477 22.378 < 2e-16 ***
## Advertising 0.1230951 0.0111237 11.066 < 2e-16 ***
## ShelveLoc.L -3.4295971 0.1082651 -31.678 < 2e-16 ***
## ShelveLoc.Q 0.3824279 0.0842896 4.537 7.62e-06 ***
## UrbanYes 0.1228864 0.1129761 1.088 0.277
## USYes -0.1840928 0.1498423 -1.229 0.220
## Population 0.0002079 0.0003705 0.561 0.575
## Income 0.0158028 0.0018451 8.565 2.58e-16 ***
## Age -0.0460452 0.0031817 -14.472 < 2e-16 ***
## Education -0.0211018 0.0197205 -1.070 0.285
## DiffPrice NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.019 on 388 degrees of freedom
## Multiple R-squared: 0.8734, Adjusted R-squared: 0.8698
## F-statistic: 243.4 on 11 and 388 DF, p-value: < 2.2e-16
Mit dem Code “Sales ~ .” können wir ein Regressionsmodell Sales gegenüber allen Variabeln erstellen.
Dank der Auswertung sehen wir, dass sieben Variabeln in diesem Modell zum Sales beitragen. Population, Education und Location (Urban und US) haben keinen Einfluss. Wir werden nun die aussagekräftigsten vertieft überprüfen.
Die Fehlermeldung bei DiffPrice kommt daher, weil DiffPrice ein abhängies Feature von anderen Feature (Differenz aus CompPrice und Price) ist, welches wir in der vorherigen Hypothese erstellt haben. Das lineare Modell erkennt, dass DiffPrice in diesem Fall keine zusätzliche Information zum Modell beiträgt. Aus diesem Grund, weist uns die Auswertung ein NA aus.
mlr.fit4 <- lm(Sales ~ Income + Advertising + ShelveLoc + Age + DiffPrice, data = Carseats)
summary(mlr.fit4)
##
## Call:
## lm(formula = Sales ~ Income + Advertising + ShelveLoc + Age +
## DiffPrice, data = Carseats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7638 -0.6955 0.0266 0.6756 3.3990
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.372761 0.222161 33.187 < 2e-16 ***
## Income 0.015903 0.001831 8.684 < 2e-16 ***
## Advertising 0.116005 0.007719 15.028 < 2e-16 ***
## ShelveLoc.L -3.415880 0.107703 -31.716 < 2e-16 ***
## ShelveLoc.Q 0.380870 0.083913 4.539 7.53e-06 ***
## Age -0.045853 0.003157 -14.522 < 2e-16 ***
## DiffPrice -0.095139 0.002660 -35.769 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.019 on 393 degrees of freedom
## Multiple R-squared: 0.8718, Adjusted R-squared: 0.8698
## F-statistic: 445.3 on 6 and 393 DF, p-value: < 2.2e-16
Durch die Auswahl der aussagekräftigsten Variabeln konnten wir ein Modell erstellen, das Sales zuverlässig erklärt. Die verwendeten Variabeln haben wir aus dem Modell entnommen, welche einen tiefen P-Wert aufweisen. In diesem Modell haben wir nur (auch in der F-Statistik) sehr tiefe p-Werte und ein ziemlich hohes R^2 von 0.87. Wir überprüfen deshalb nun die Residuen.
par(mfrow = c(2,2))
plot(mlr.fit4)
Im Q-Q Plot wird bestätigt, dass es sich um ein lineares Modell handelt. Der erste Plot und das folgende Histogram bestätigen ebenfalls, dass die Residuen normalverteilt und einen Mittelwert von 0 aufweisen.
hist(mlr.fit4$residuals,
main = "Verteilung der Residuen", col = "lightblue")
Hypothese: Sales lässt sich durch die Verwendung von mehreren Features zuverlässiger vorhersagen.
Unser letztes Modell hat gezeigt, dass mit einer sorgfältigen Auswahl der Variabeln, Sales gut und zuverlässig erklärt werden kann. In der Praxis könnten noch mehr Kombinationen (Variabeln, Interaktionen, exponentielle Zusammenhänge etc…) untersucht werden. In dieser Abgabe lassen wir das weg, weil es sich um Wiederholungen der gezeigten Funktionen handelt.
Thema: Logistische Regression
Um logistische Regressionen durchzuführen, benötigen wir binäre Variabeln. Im ersten Schritt ersetzen wir Yes/No von Urban durch 1/0.
Carseats <- Carseats %>%
mutate(Urban_binaer = recode(Urban, "No" = 0, "Yes" = 1))
ggplot(Carseats, aes(x=Sales, y=Urban)) +
geom_point(color = "blue") +
labs(x = "Sales",
y = "Urban",
title = "Sales vs Urban",
subtitle = "Carseats Datensatz")
Da wir in den Plots zur logistischen Regression die Sigmoid Funktion einzeichnen möchten, verwenden wir diese Darstellung bei sämtliche logistischen Regression.
Dank der visuellen Analyse von Sales und Urban wird ersichtlich, dass Urban:Yes in beide Richtungen stärker streut als Urban:No. Entsprechend wird das Modell kaum verlässliche Werte ergeben.
logreg.fit1 <- glm(Urban_binaer ~ Sales, data = Carseats , family = binomial)
summary(logreg.fit1)
##
## Call:
## glm(formula = Urban_binaer ~ Sales, family = binomial, data = Carseats)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6032 -1.5464 0.8288 0.8396 0.8737
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.96122 0.31250 3.076 0.0021 **
## Sales -0.01197 0.03883 -0.308 0.7578
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 485.25 on 399 degrees of freedom
## Residual deviance: 485.16 on 398 degrees of freedom
## AIC: 489.16
##
## Number of Fisher Scoring iterations: 4
Wie wir aufgrund der Visualisierung erkannt haben, ergibt ein logistisches Regressionsmodell einen grossen P>|z| Wert. Somit gibt es keinen logistischen Zusammenhang zwischen Sales und Urban.
Gegenüber der lineare Regression haben wir den Parameter “family = binomial” erfasst.
Um ein Beispiel für eine logistische Regression zeigen zu können, untersuchen wir auch Sales und ShelveLoc
Zuerst für alles ShelveLoc, danach haben die Location “Medium” entfernt.
Carseats %>%
ggplot(aes(x = Sales, y = ShelveLoc)) +
geom_point(color = "blue") +
labs(x = "Sales",
y = "ShelveLoc",
title = "Sales vs ShelveLoc",
subtitle = "Carseats Datensatz")
Carseats_ex_medium <- Carseats %>% filter(ShelveLoc != "Medium")
Carseats_ex_medium %>%
ggplot(aes(x = Sales, y = ShelveLoc)) +
geom_point(color = "blue") +
labs(x = "Sales",
y = "ShelveLoc",
title = "Sales vs ShelveLoc",
subtitle = "Carseats Datensatz")
Im ersten Plot wurden drei Variabeln dargestellt. Für ein gutes Vorzeigebeispiel verwenden wir nur Good und Bad.
Carseats_ex_medium <- Carseats_ex_medium %>%
mutate(ShelveLoc_binaer = recode(ShelveLoc, "Bad" = 0, "Good" = 1))
head(Carseats_ex_medium)
## Sales Price CompPrice Advertising ShelveLoc Urban US Population Income Age
## 1 9.50 120 138 11 Bad Yes Yes 276 73 42
## 2 11.22 83 111 16 Good Yes Yes 260 48 65
## 5 4.15 128 141 3 Bad Yes No 340 64 38
## 6 10.81 72 124 13 Bad No Yes 501 113 78
## 8 11.85 120 136 15 Good Yes Yes 425 81 67
## 11 9.01 100 121 9 Bad No Yes 150 78 26
## Education DiffPrice Urban_binaer ShelveLoc_binaer
## 1 17 -18 1 0
## 2 10 -28 1 1
## 5 13 -13 1 0
## 6 16 -52 0 0
## 8 10 -16 1 1
## 11 10 -21 0 0
Wir haben “Medium” weggelassen und wieder eine binäre Klassifizierung erstellt.
logreg.fit2 <- glm(ShelveLoc_binaer ~ Sales, data = Carseats_ex_medium , family = binomial)
summary(logreg.fit2)
##
## Call:
## glm(formula = ShelveLoc_binaer ~ Sales, family = binomial, data = Carseats_ex_medium)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3930 -0.5021 -0.1506 0.5015 2.5892
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.0248 0.8714 -6.914 4.71e-12 ***
## Sales 0.7566 0.1082 6.990 2.76e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 250.25 on 180 degrees of freedom
## Residual deviance: 136.49 on 179 degrees of freedom
## AIC: 140.49
##
## Number of Fisher Scoring iterations: 5
Ein tiefer P-Wert weist auf eine hohe Zuverlässigkeit vom Modell hin.
Im folgenden Plot wird ersichtlich, dass das Modell bei allen Läden mit Sales > 7.96 den ShelveLoc als Good klassifiziert. Die Berechnung des Wertes haben wir im folgenden Code Chunk auskommentiert.
ggplot(Carseats_ex_medium, aes(x=Sales, y=ShelveLoc_binaer)) +
geom_point(color = "blue") +
labs(x = "Sales",
y = "ShelveLoc",
title = "Sales vs ShelveLoc",
subtitle = "Carseats Datensatz") +
stat_smooth(method="glm", formula = y ~ x, color="red",
se = FALSE, method.args = list(family=binomial)) +
geom_vline(xintercept = 7.9633)
# Berechnung der Decision Boundary
#-logreg.fit2$coefficients[1] / logreg.fit2$coefficients[2]
Um die Metriken des Modelles zu berechnen, erstellen wir eine Confusion Matrix.
# zuerst erstellen wir eine Prediction, mit den Testdaten
# Predicte ShelveLoc_binaer Werte
logreg.fit2.probs <- predict(logreg.fit2, type = "response")
# Erstellt einen Vektor mit 181 Null Werten
logreg.fit2.pred <- rep(0, 181)
# Uebreschreibt alle Werte die groesser sind als 0.5 als 1
logreg.fit2.pred[logreg.fit2.probs > .5] = 1
# Erstellen eine Kreuztabelle (Confusions Matrix)
table("Actual Class" = Carseats_ex_medium$ShelveLoc_binaer, "Predicted Class" = logreg.fit2.pred)
## Predicted Class
## Actual Class 0 1
## 0 82 14
## 1 19 66
Die Spalten entsprechen den tatsächlichen Werten und in den Zeilen befinden sich die vorhergesagten Werte. Lesebeispiel: Das Modell hat 19 Werte als 0/Bad klassifiziert, die aber tatsächlich ein 1/Good gewesen sind.
Die Accuracy entspricht dem Anteil korrekter Prognosen an allen Prognosen.
# Accuracy Berechnen
acc <- (82 + 66) / (82 + 19 + 14 + 66)
acc
## [1] 0.8176796
Precision = True Positiv / (True Positiv + False Positiv) Preicision ist, Anteil richtig positiv an allen Positiv
Precision entspricht der Anteil korrekter positiven/1 Prognosen an allen positiven Prognosen.
# Precision Berechnen
precision <- 66 / (66 + 14)
precision
## [1] 0.825
Recall = True Positiv / (True Positiv + False Negativ) Recall ist, Anteil richtig positiv an allen richtigen Vorhersagen
Recall / Sensitivity entspricht der Anteil korrekter positiven/1 Prognosen an allen tatsächlich positiven Werten.
# Recall berechnen
recall <- 66 / (66 + 19)
recall
## [1] 0.7764706
Mit unseren verwendeten Modellen konnten wir aufzeigen, dass Sales nicht von einzelnen Variabeln, aber von mehreren zusammenhängt. Dank der multiplen linearen Regression konnten wir zuverlässige Einflussfaktoren aufzeigen.
Dank der logistischen Regression konnten wir darstellen, dass Urban keinen Einfluss hat. Mit dem Weglassen von Medium konnten wir für ShelveLoc ein zuverlässiges Modell erstellen.